library(dplyr)
# For Model fitting
library(lme4)
library(purrr)
# For diagnostics
library(performance)
# For adding new columns
library(tibble)
# For highlighting min value in a table
library(reactablefmtr)
# For boxcox function
library(caret)# Load data
sys.source("./codes/scripts/code_join_data_full_dataset.R", envir = knitr::knit_global())# Load functions
## Models
sys.source("./codes/functions/functions_models.R", envir = knitr::knit_global())
## Inference
sys.source("./codes/functions/function_anova_and_tukey_tables.R", envir = knitr::knit_global())
## Plot
sys.source("./codes/functions/function_masomenos_plots.R", envir = knitr::knit_global())Select performance and traits variables
data_for_models <-
data_complete %>%
# Calculate nitrogen use efficiency column
# I followed Leaf traits explaining the growth of tree
# species planted in a Central Amazonian disturbed area
mutate(photo_nitrogen_use = (amax*sla_cm2_g*0.1)/perc_n) %>%
# select variables that are going to be used in the models
select(id, spcode, treatment,nfixer, init_height,
# Plant preformance
total_biomass, above_biomass, below_biomass, agr,
# Traits
amax, gs, wue, d13c, d15n, photo_nitrogen_use) %>%
mutate(id = factor(id))
data_for_models$nfixer <- relevel(data_for_models$nfixer,ref = "fixer" )
levels(data_for_models$nfixer)[1][1] "fixer"
Box-Cox transformation
# Nonnormality and heterogeneity of variance were corrected by applying box-cox transformation
box_cox <-
data_for_models %>%
# Transforme only the traits
select(5,6:9,10:15) %>%
purrr::map(., ~ BoxCoxTrans(.x))Add variables transformed to the dataset
data_for_models_transformed <-
cbind(data_for_models,
# Plant's performance
total_biomass_transformed = predict(box_cox$total_biomass, data_for_models$total_biomass),
above_biomass_transformed = predict(box_cox$above_biomass, data_for_models$above_biomass),
below_bioomass_transformed = predict(box_cox$below_biomass, data_for_models$below_biomass),
agr_transformed = predict(box_cox$agr, data_for_models$agr),
# Traits
amax_transformed = predict(box_cox$amax, data_for_models$amax),
gs_transformed = predict(box_cox$gs, data_for_models$gs),
wue_transformed = predict(box_cox$wue, data_for_models$wue),
# d13 and d15 where not transformed because the data has negative values
photo_nitrogen_use_transformed = predict(box_cox$photo_nitrogen_use, data_for_models$photo_nitrogen_use),
# Covariate
init_height = predict(box_cox$init_height, data_for_models$init_height)) %>%
# Remove original variables (non-transformed)
select(id, spcode, treatment, nfixer, init_height, everything(),
-c(5,6:9,10:12,15)) Models: Questions 1 and 2
\[response\sim treatment*fixer\ + initial\ height + random( 1|\ specie)\]
# Take response variables' names
response_vars_q1_q2 <- set_names(names(data_for_models_transformed)[6:(ncol(data_for_models_transformed))])model_list_q1_q2 <- map(response_vars_q1_q2, ~ mixed_model_1(response = .x, data = data_for_models_transformed))Models: Questions 3
\[performance\sim treatment*fixer*trait\ + initial\ height + random( 1|\ specie)\]
# This funtion takes two columns a create a model formula with all the possible
# combinations
model_combinations_formulas <- function(y_var, x_var){
variables <- crossing(y_var, x_var)
pattern <- c('\\+.*|nfixer|treatment|[[:punct:]]| ')
# Model
models <- paste0(variables$y_var,"~treatment*nfixer*",
variables$x_var,"+init_height+(1|spcode)")
formulas <- map(models, as.formula)
# Add name to the list
names(formulas) <- stringr::str_replace_all(formulas, pattern, replacement = '')
return(formulas)
}Scale preditors
data_for_models_transformed_scaled <-
data_for_models_transformed %>%
# scale pedictors
mutate(across(c(5, 10:15), scale))# Select traits (x_vars)
traits_names <- set_names(names(data_for_models_transformed_scaled)
[c(6,7, 12:15)])
# Select plants performance vars (y_vars)
performance_names <- set_names(names(data_for_models_transformed_scaled)[8:11])
models_lmer_formulas <- model_combinations_formulas(performance_names, traits_names)model_list_q3 <- map(models_lmer_formulas,
~ lmer(.x, data = data_for_models_transformed_scaled))Validation plots performance package
Models for questions 1,2
map(model_list_q1_q2, check_model)$d13c
$d15n
$total_biomass_transformed
$above_biomass_transformed
$below_bioomass_transformed
$agr_transformed
$amax_transformed
$gs_transformed
$wue_transformed
$photo_nitrogen_use_transformed
Models for question 3
map(model_list_q3, check_model)$`abovebiomasstransformed~amaxtransformed`
$`abovebiomasstransformed~d13c`
$`abovebiomasstransformed~d15n`
$`abovebiomasstransformed~gstransformed`
$`abovebiomasstransformed~photonitrogenusetransformed`
$`abovebiomasstransformed~wuetransformed`
$`agrtransformed~amaxtransformed`
$`agrtransformed~d13c`
$`agrtransformed~d15n`
$`agrtransformed~gstransformed`
$`agrtransformed~photonitrogenusetransformed`
$`agrtransformed~wuetransformed`
$`belowbioomasstransformed~amaxtransformed`
$`belowbioomasstransformed~d13c`
$`belowbioomasstransformed~d15n`
$`belowbioomasstransformed~gstransformed`
$`belowbioomasstransformed~photonitrogenusetransformed`
$`belowbioomasstransformed~wuetransformed`
$`totalbiomasstransformed~amaxtransformed`
$`totalbiomasstransformed~d13c`
$`totalbiomasstransformed~d15n`
$`totalbiomasstransformed~gstransformed`
$`totalbiomasstransformed~photonitrogenusetransformed`
$`totalbiomasstransformed~wuetransformed`
Save lists with the models
saveRDS(model_list_q1_q2, file = "./processed_data/models_q1_q2.RData")
saveRDS(model_list_q3, file = "./processed_data/models_q3_3_way_interaction.RData")